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ABSTRACT 

We review recent progress in understanding the role of chaos in influencing the structure 
and evolution of galaxies. The orbits of stars in galaxies are generically chaotic: the 
chaotic behavior arises in part from the intrinsically grainy nature of a potential that is 
composed of point masses. Even if the potential is assumed to be smooth, however, much 
of the phase space of non-axisymmetric galaxies is chaotic due to the presence of central 
density cusps or black holes. The chaotic nature of orbits implies that perturbations will 
grow exponentially and this in turn is expected to result in a diffusion in phase space. 
We show that the degree of orbital evolution is not well predicted by the growth rate 
of infinitesimal perturbations, i.e. by the Liapunov exponent. A more useful criterion is 
whether perturbations continue to grow exponentially until their scale is of order the 
size of the system. We illustrate these ideas in a potential consisting of N fixed point 
masses. Liapunov exponents are large for all values of iV, but orbits become increasingly 
regular in their behavior as N increases; the reason is that the exponential divergence 
saturates at smaller and smaller distances as N is increased. The objects which lend 
phase space its structure and impede diffusion are the invariant tori. In the triaxial 
potentials we discuss, a large fraction of the tori correspond to resonant (thin) orbits 
and their associated families of regular orbits. These tori are destroyed by perturbations 
to the potential. When only a few stable resonances remain, we find that the phase space 
distribution of an ensemble of chaotic orbits evolves rapidly toward a nearly stationary 
state. This mixing process is shown to occur on timescales of a few crossing times in 
triaxial potentials containing massive central singularities, consistent with the rapid 
evolution observed in A^-body simulations of galaxies with central black holes. 



1 INTRODUCTION 



The role of microscopic chaos in producing macroscopic re- 
laxatiom of dynamical systems has received considerable at: 
tention In Lliti held of sLaLlsLlcal meclianlcs (LebowlL^ 1993 



on which the discrete component of the force is important 
is much longer than the age of the Universe. At the same 
time , it has been shown (Miller 1964, Gurzadyan & Savvidy 



Gaspard 1998). However, the role of chaos in the relaxation 



of stellar systems is relatively less well understood. The last 
five years have witnessed the development of a number of 
new techniques for probing the complexities of phase space 
in realistic galactic potentials. This work has led to a greater 
understanding of the importance of stochastic orbits and 
their role in driving dynamical evolution. 

The gravitational forces on a star in a galaxy can be 
broken up into two components: a rapidly varying compo- 
nent that arises from the discrete distribution of stars, and 



1986) that trajectories in the A^-body problem are generi- 
cally chaotic, with rates of exponential divergence that ap- 
pear to remain large even for large A'^. Thus, according to at 
least one definition of chaos, the orbits of stars in galaxies 
should never tend toward the regular motion expected in 
smooth potentials. 

The exponential divergence of trajectories gives rise to 
a diffusion. The approach of a non-stationary distribution of 
phase points to a stationary one via chaotic di ffusion is re- 
ferred to as "chaotic mixing" or just "mixing" (Kandrup & 
Mahon 1994). In the context of non-equilibrium statistical 



a smoothly varying component that arises from the overall 
mass d istribution. The importance of the discrete compo- 
nent of the force relative to the smooth component is usually 



assumed proportional to ~ Win N/N, the ratio of dy nami- 
cal to two-body relaxation times ( Ghandrasekhar 1943). The 



implication is that, for galaxies (TV ~ 10^^), the time scale 



mechanics, mixing to an invariant distribution is regarded 
as a legitimate relaxation process ( ^ebowitz 199S , Gaspard 
199S). However, in the context of gravitational systems, the 
issue of whether or not there is a connection between the 
time scale for exponential divergence of adjacent trajecto- 
ries and the time scale for macroscopic relaxation has been 
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a matter of much debate ([Gurzadyan fc Savvidy 1986, |Kan- 
drup idM [Kandrup fc Smith 1994 pcggie 1991| , poodm an 
et al. l' )9%. 



The question of how and under what conditions chaos 
will result in evolution of a stellar system is the subject of 
this article. In §2 we present two idealized models: a set of 
A'^ fixed mass points, representing stars in a galaxy; and a 
smooth potential with a single central point mass, repre- 
senting a galaxy with a nuclear black hole. Motion in both 
potentials is generically chaotic, with similar values of the 
Liapunov exponent. We show that the Liapunov exponent 
- which measures only the infinitesimal growth of pertur- 
bations - is not always useful for predicting the degree of 
macroscopic evolution. A more useful criterion is whether 
perturbations continue to grow exponentially until their size 
is of order the size of the system. The kinds of structures that 
can impede diffusion of stochastic orbits arc discussed in §3. 



harmonic oscillator, and one expects to see corresponding 
changes in the behavior of orbits. 

The evolution of orbits in Model 1 is illustrated in Fig. 
1, as a function of the number of particles making up the 
potential. For each A'^, a set of 32 random realizations of the 
particle positions were generated; a single orbit was then 
integrated in each of the 32 corresponding potentials. The 
evolution of the orbits were measured in several ways. The 
mean Liapunov exponent a (averaged over the 32 ensem- 
bles) was found to increase with A'^ until TV ~ 100 then level 
off, at a value aTcr ~ 5, with Tc- defined as one-half of the 
long-axis orbital period in the smooth potential. Thus these 
trajectories are locally unstable on a remarkably short time 
scale, just a fraction of a crossing time, and this characteris- 
tic time shows no tendency to decrease with increasing A*' for 
A' as large as ~ 10^ . This result is similar to the well-known 



exponential instability of the full A^-body problem (Miller 



In §4 W 3 simulate the approach to an invariant distribution 1964 Gurzadyan fc Savvidy 1986 , Kandrup 1990 , Kandrup 



of start in galactic potentials; we find that evolution to a & Smith 1991, [Heggie 1991 , Goodman et al. 1993) 



stationary state can take place in little more than a crossing 
time if the phase space is globally chaotic. Such evolution is 
consistent with the rapid, self-consistent evolution observed 
in numerical simulations of galaxies containing nuclear black 
holes. 



2 THE LIAPUNOV EXPONENT AND OTHER 
MEASURES OF CHAOS 

Exponential divergence of nearby trajectories is a common 
property of dynamical systems. This divergence is most com- 
monly expressed in terms of the Liapunov exponent, which 
measures the mean e-folding rate of an infinitesimal per- 
turbation averaged over an infinite time interval. Because 
the perturbation is assumed to remain small, Liapunov ex- 
ponents measure only the rate of divergence in the imme- 
diate vicinity of the unperturbed orbit; they do not nec- 
essarily contain any information about the non-linear, or 
macroscopic, evolution. One might nevertheless expect the 
magnitude of the Liapunov exponent to predict, in some ap- 
proximate sense, the degree to which chaos will induce finite 
changes in the structure of an orbit after a fixed time. This 
turns out not always to be the case, as we now illustrate. 

Figures 1 and 2 summarize the results of test-particle 
integrations in two time-independent potentials. Model 1 
(Fig. 1) consists of A'^ point masses m distributed randomly 
and uniformly within a triaxial ellipsoid; the total mass 
M = Nm remains fixed (M — 1) as TV is varied, as do 
the axis lengths {a = l,b — 0.75, c = 0.5) of the ellipsoid. 
Model 2 (Fig. 2) is the smooth representation of Model 1, 
i.e. a homogeneous ellipsoid, to which has been added a sin- 
gle central point of mass Mh- The parameter TV of Model 1 
may be interpreted as the number of stars that make up a 
galaxy of fixed mass, and the parameter Mh of Model 2 as 
the mass of a central black hole, expressed in units of the to- 
tal galaxy mass. The quantities 1/TV and Mh play the role of 
perturbation parameters; as they are increased, the poten- 
tial departs more and more from that of the integrable, 3D 



As a second measure of the orbital evolution, the RMS 
variations in the x, y and z "actions" of the orbits were 
computed over ~ 20 periods; these "actions" were defined 
in terms of the frequencies of motion in the smooth poten- 
tial, e.g. Ja: = Ex/oJx = (f^ -\- LOxX^) / {"^.tOx), and would be 
precisely conserved in the limit of zero perturbation. Figure 
1 shows that there is a smooth decrease in the amplitude of 
5,Jx as TV is increased - in other words, the orbits approach 
more and more closely, in their macroscopic behavior, to 
that of the integrable orbits even though they remain locally 
unstable to a degree (as measured by cr) that is nearly inde- 
pendent of TV. Plots of the trajectories (Fig. 1) confirm that 
the orbits for large TV are similar to the Lissajous figures 
expected in the 3D harmonic oscillator. 

This apparent paradox - large cr, but nearly regular mo- 
tion - is reconciled in the lower right panel of Fig. 1, which 
shows the variation with time of the finite separation be- 
tween two initially nearby trajectories, for four values of TV. 
The divergence is initially exponential in all cases, with a 
characteristic time ~ 0.2Tcr, consistent with the measured 
values of a. However the exponential divergence eventually 
saturates, after which the separation slows. This saturation 
has no effect on the Liapunov exponents which are calcu- 
lated assuming that the separation remains infinitesimally 
small. Furthermore, the saturation occurs sooner for larger 
TV. When TV < 10^, the exponential divergence continues un- 
til the separation is of order the size of the system, while for 
TV ~ 10^, saturation occurs at a separation of only ^ 10~^, 
much smaller than the system size. 

How can the separation between two, locally unstable 
orbits saturate at a small amplitude, given that neither orbit 
"knows" about the other orbit? The answer must be that 
both orbits are confined to the same, restricted region of 
phase space; saturation occurs when the separation between 
them is of order the width of this region. Further systematic 
increase in their separation would require that one of the 
trajectories "break out" into another such region, and such 
events (for large TV) apparently occur at a much lower rate 
than the divergence described by the Liapunov exponent. 
The fact that the exponential divergence saturates sooner 
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Figure 1. Evolution of orbits in a potential consisting of A'^ fixed point masses, m = 1/A'^, djstributed randomly and uniformly in an 
ellipsoidal volume. For each TV, 32 random realizations of the particle positions were generated and a single orbit was integrated in each 
of the corresponding potentials. The upper left panel shows mean values of the Liapunov exponent a for each ensemble, and the RMS 
variation in the "action" Jx over ~ 20 orbital periods (see text), a remains large while 5Jx falls off with A^; thus the behavior of a typical 
orbit becomes more and more regular as TV increases. This is verified in the lower left panel which shows typical orbits for four values of 
TV. The lower right panel displays the finite separation between initially adjacent trajectories. While the divergence takes place at the 
rate predicted by the Liapunov exponents at early times, for large TV the separation saturates at a value much smaller than the size of 
the system. 



for larger TV suggests that the width of the confining regions 
decreases with increasing N - although our experiments do 
not allow us to estimate the precise TV-dependence. The ap- 
parently unconfined evolution for TV < 10^ suggests that 
phase space for such small TV is "globally chaotic", with es- 
sentially no confining barriers. 

Figure 1 suggests the way in which the equations of 
motion in an TV-body potential tend toward those of the 
corresponding smooth potential as TV increases. Although 
some measures of chaos - e.g. the Liapunov exponents - 
remain large even for large TV, others - e.g. the change in 
the actions - tend to zero. 

These results present an interesting contrast with those 
from Model 2 (Fig. 2) . Here a single orbit was integrated for 
each value of Alh, the mass of the central "black hole," in 
the otherwise smooth, harmonic-oscillator potential. Both a 



and SJx were now found to vary systematically with pertur- 
bation parameter Mh ; thus the orbits tend toward regularity 
with decreasing Mh as measured both by their infinitesimal 
and finite-amplitude behaviors. But the lower right-hand 
panel of Fig. 2 reveals essential similarity with the behav- 
ior of the orbits in Model 1. For large black hole masses, 
Mh > 0.02, divergence continues at roughly the Liapunov 
rate until the orbits are separated by of order the size of 
the system; while for Mh < 0.01, the evolution of the sepa- 
ration is more complex, with periods of stagnation followed 
by sudden jumps (when the trajectory passes sufficiently 
close to the central point). Once again, we postulate that 
the phase space of potentials with A'lh > 0.01 is globally 
stochastic, with few impedime nts to the motion - a co njec- 
ture that was actually verified (Valluri & Merritt 199J) in a 



different class of triaxial potentials (see §3). For smaller Mh, 
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Figure 2. Like Fig. 1, for orbits in tlie potential of a uniform triaxial ellipsoid of unit mass ■^ith an added central point of mass M^. 
A single orbit was integrated for each value of Mj, and its properties recorded. Both the Liapunov exponent a and the variation in the 
"action" SJx fall off with decreasing M^. However, as in Fig. 1, the trajectories can mimic regular orbits for times much longer than cr~^ 
when M/i is small. The lower right panel again shows why: for small Af^, the evolution of the finite separation between initially adjacent 
trajectories is not well predicted by the Liapunov exponent. 



the orbital evolution consists of a sequence of nearly-random 
t ransitions from one, app roximately regular orbit to another 
( Gerhard & Binney 1985 ) ; these transitions become rarer as 
Mh is reduced. 

These experiments suggest that Liapunov exponents are 
not good predictors of the macroscopic behavior of orbits. 
Even trajectories with large a can exhibit nearly regular 
behavior over time scales much longer than . In decid- 
ing whether stochasticity is likely to be important for the 
behavior of orbits, the crucial question is not the value of 
o", but whether the exponential divergence continues until 
separations of order the size of the system are reached. If 
the separation saturates at some value much smaller than 
the system size, orbital stochasticity may not be of much 
consequence even if a is large. 

The sorts of structures that can impede the diffusion 
of stochastic orbits are described in more detail in the next 
section. Here we note that even some very simple dynamical 



models can exhibit the basic properties that we have de- 
scribed. For instance, the motion of particles in a "Lorentz 
gas," a fixed 2D array of cylindrical scatters, is generically 
stochastic, with Liapunov times of order the time between 
collisions. However in the close-packed limit, trajectories are 
blocked by the finite size of the scatterers, causing them to 
r emain confine d to narrow regions for long periods of time 
( [Gaspard 199^ ). Thus the exponential instability has little 
consequence for the macroscopic character of the motion. 



3 PHASE SPACE STRUCTURE OF TRIAXIAL 
GALAXIES 

In an integrable potential with N degrees of freedom (DOF), 
all the trajectories have exactly isolating integrals and are 
confined to A^'-dimensional tori in phase space. Such a torus 
is defined by the A*' actions J that fix its cross-sectional ar- 
eas. Motion around this torus occurs at a rate determined by 
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a constant frequency vector {uji, u!2, (^n)- Realistic poten- 
tials with more than 1 DOF are rarely integrable and the 
motion is more complex. While the KAM theorem guar- 
antees that most of the original tori will persist when an 
integrable system is slightly perturbed, even under small 
perturbations a large part of the phase space will be influ- 
enced by resonances. A resonant torus is one which satisfies 
(one or more) resonant conditions of the form n.Lj= be- 
tween the A'' fundamental frequencies. Such tori are dense 
in the phase space of the integrable potential but typically 
only the lowest order resonances have a signi ficant influence 
on the motion i n the perturbed potential ( Lichtenberg & 
Lieberrian f992). In the vicinity of a stable resonant torus, 
motion is still regular, but the orbits have shapes determined 
by the order of the resonance - often very different from the 
shapes of orbits in the integrable potential. In the vicinity 
of unstable resonant tori, trajectories are usually chaotic. 

In a system with 2 DOF, the resonant tori are closed 
periodic orbits defined by a single frequency wo = n2UJi = 
niLJ2- In 3 DOF the resonant tori obey a condition on the 
three fundamental frequencies of the form niui + n2UJ2 + 
naus = 0. Such a relation does not imply that the orbit is 
closed, as in 2 DOF, but rather th at it is thin, densely fil l- 
ing a sheet in configuration space (Vlerritt & Valluri 1999). 



The resonance condition can be used to reduce the number 
of independent frequencies by one; the two remaining fre- 
quencies then describe the rate of rotation around the 2D 
reduced torus ( [Carpintero fc Aguilar 199^ , Merritt & Val- 
luri 199^. In 3 DOF systems as in 2 DOF systems, resonant 
tori are the regions where perturbation expansions fail, and 
where the the global structure of phase space is expected 
to be changed. Thus the thin orbits are expected to play a 
similar role in 3 dimensions to the role played by periodic 
orbits in 2 dimensions. 

Families of thin orbits may exist even in an integrable 
potential if it co ntains "primary resonances" th at divide up 
the phase space ( Lichtenberg fc Lieberman 1992|). This is the 



case in the w ell-known "perfect ellipsoid" ( |Kuzmin 1973 



Zeeuw L985): the 1 : 1 closed orbits in the principal planes 
generate families of thin tube orbits. In non-integrable po- 
tentials, one expects to find thin orbits and their associated 
families throughout phase space. 

One of the most useful tools to be employed recently 
in the study of galactic potentials is the fre quency analy- 



sis technique (NAFF) develop ed by Laskar (Laskar 1990 
Laskar et al 



1992 



Laskar 1998). Laskar's technique exploits 



the fact that regular orbits are quasiperiodic: Cartesian co- 
ordinates like x{t) and Vx{t) can be expressed as Fourier 
series in terms of the three fundamental frequencies on the 
torus. While this basic principle has been used to study or- 
bits in the past (Binney & Spergel 1984), Laskar showed 



that the accuracy with which the individual frequency com- 
ponents can be extracted is greatly improved by employing a 
Fourier filtering function and by orthogonalizing successive 
frequency components. Once the entire frequency spectrum 
is obtained, the three fundamental frequencies (tJi, tJ2, 'i's), 
which appear as linear combinations in each line, may be 



identifi ed using an integer programming algorithm (Valluri 
Men itt 1998 ) . All the lines in the spectrum are then im- 



mediately expressible in terms of the fundamental frequen- 
cies, giving a map between the actio n-angle and Cartesian 
coordinates ( Valluri fc Merritt 1999b ). 

Strictly speaking, only regular orbits are restricted to 
tori and amenable to Laskar's technique. Nonetheless, as 
discussed above, stochastic trajectories may mimic regular 
orbits for long periods of time. On short time scales, such 
an orbit has a frequency spectrum which mimics a quasi- 
periodic series. As the orbit diffuses through phase space , 
its frequency spectrum will change. Laskar ( Laskar 199^ ) 
showed that over a fixed interval of time (AT) the change 
in the frequency of the leading term in the spectrum, Auj = 
tj(T)— aj(r-|-Ar)|, is a good measure of the rate of diffusion 
of a chaotic orbit in phase space. Note that Laskar's Aoj, 
unlike the Liapunov exponent, measures a finite excursion 
in phase space. In this sense it is similar to the 5 J parameter 
defined above. 

Since Laskar's technique maps the structure of phase 
space in the frequency domain, it is ideally suited to iden- 
tifying the resonant t ori. In 2 DOF galactic potentials, Pa- 
paphillipou fc Laskar (1996) showed that in a map of uj\/uo2 
versus a third parameter, resonant families appear as regions 
where u)-i/u)2 = const over some set of orbits. In 3DOF sys- 
tems, a "frequency map" m ay be obtained by plotting th e 
ratios ui/uz versus uJ2/(jJz ( Papaphillipou fc Laskar 1998). 



Here the resonant families show up as lines of constant slope. 
Thus by constructing a frequency map of an ensemble of or- 
bits chosen from a regular grid of starting values, it is pos- 
sible to locate the resonances which significantly affect the 
structure of phase space. 

Another useful device is t he "diffusion map" : a plot of 
the chaotic diffusion rates Aoj (Laskar 1993). We calculated 



diffusion maps for a family of galaxy models which are tri- 
axial ge neraliz ations of the spherical models described by 
Dehnen (1993) and others. These models provide a good fit 



to the inner part of the light distributions of real galaxies 
and have a mass density given by 

(3-7)M 



p(m) 



with 



4mabc 



1 



-(4-7) 



< 7 < 3 



y_ 

b2 



+ 



a > fe > c > 0, 



(1) 



(2) 



and M the total mass 
with axis ratios a : b 



The mass is stratified on ellipsoids 
c; the X and z axes are the long 
and short axes respectively. The parameter 7 determines 
the slope of the central, power-law density cusp. 

Figures 3a, b show the diffusion maps at one energy in a 
triaxial Dehnen potential with 7 = 0.5, with and without a 
central black hole. Orbits were integrated starting at rest on 
the equipotential surface just inside the half-mass radius of 
the model. The diffusion rate (Ali;) was plotted in grey scale 
on the equipotential surface. Rapidly diffusing (chaotic) or- 
bits are coded in black while regular orbits (|Aa;| ~ 0) are 
shaded white. The intensity of the grey scale is proportional 
to the logarithm of the stochastic diffusion rate as measured 
by |Ati;|. The resonant tori and their associated families ap- 
pear as narrow white bands and are marked in the figure 
by their defining integers n\,n2,nz. The white bands are 
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Figure 3. Resonances in triaxial potentials. The mass model in (a) has a weak (7 = 0.5) cusp and no black hole; in (b) the black hole 
contains 0.3% of the total mass. Both equipotential surfaces lie close to the half-mass radius. The grey scale measures the degree of 
stochasticity of orbits started with zero velocity on the equipotential surface. Stable resonance zones - the white bands in (a) and (b) 
— are labeled by the order {ni,n2,n3) of the resonance. Panels (c) and (d) show the pericenter distance A of a set of 10'^ orbits with 
starting points lying along the heavy solid lines in (a) and (b). 



flanked by thin dark ones which are the stochastic layers 
marking the transitions between resonant families. Some of 
these resonances - the (2, -1) banana; the (4, -3,0) pretzel; 
and the (3,0,-2) fish - are 2D resonances that correspond 
to periodic orbits in one of th e principal planes. Miralda- 
Escude & Schwarzschild (1989) call such orbits "boxlets". 
However most of the resonant families cannot be associated 
with a periodic orbit. 

A striking feature of Fig. 3a is the large number of nar- 
row but distinct resonance zones. The reason for the nar- 
rowness is suggested by Fig. 3c which shows the distance 
of closest approach to the potential center of a set of orbits 
whose initial conditions lie along the heavy curve in Fig. 3a. 
As one passes through a stable resonance, the orbital peri- 
center reaches a maximum on the resonance where the orbit 
has zero thickness. The width of the resonance band is set by 



how thick an orbit can become before it is rendered stochas- 
tic - this happens when the distance of closest approach 
to the destabilizing center is sufficiently small. In potentials 
with a steep central force, this destabilization is attributable 
to gravitational deflections which occur when the trajectory 
passes near the center. In the case of the 7 — 0.5 model of 
Fig. 3a, the force is analytic all the way to r = 0. Here, the 
orbits become chaotic because they come arbitrarily close to 
an unstable "separatrix" layer which marks the transition 
between one resonant family and the next. Once a family 
is wide enough to pass through the center it can no longer 
maintain the same fundamental "shape" properties (i.e. it 
can no longer be characterized by the same (711,712,^3) val- 
ues). This results in a transition to a new family of orbits. 
This transition from one resonant family to another is rem- 
iniscent of the transition from tube to box orbits, which 
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occurs even in integrable potentials ( |de Zeeuw 1985| ). The 
lowest-order resonant orbits, have simple shapes and the 
largest central "holes" allowing a large family of associated 
regular orbits. Higher-order resonances have more complex 
shapes and pass nearer the center; their associated families, 
and their resonance bands in Fig. 3a, are correspondingly 
smaller. 



"imperfect ellipsoidal" potential (Merritt & Valluri 199(:) 



Fif ure 3b shows a diffusion map for a 7 = 0.5 model 



with a central black hole of mass ~ 0.003 A/. In this 
figure the stochastic regions between the resonant families 
are significantly wider. Figure 3d shows that the black hole 
destabilizes orbits well before they become thick enough to 
pass through the center - somewhat dif ferent from the stan- 
dard picture ( [Gerhard fc Binney 1985 ) in which box orbits 
are destabilized by passing arbitrarily close to a central black 
hole. 

As the mass of a central point is increased, more and 
more of the resonant tori are rendered unstable, and the 
width of the stochastic layers and their degree of overlap ar e 
increased. Figure rotation has a similar effect (Valluri 1999). 



One finds (V^alluri fc Merritt 1998, Merritt fc Valluri 199£ ) 



that the phase space of box-like orbits (orbits with station- 
ary points) becomes "globally chaotic" when such perturba- 
tions are sufficiently large - there are essentially no stable 
invariant tori left. For instance, in non-rotating triaxial po- 
tentials, global chaos ensues when the central point contains 



~ 10"'^ of the galaxy mass within the equipotential surface. 
In the globally chaotic regime, there are few barriers to the 
motion and one expects orbital evolution to be very rapid 
and extensive. Figure 2 confirms this prediction in the case 
of motion in a homogeneous ellipsoid. 



CHAOTIC MIXING AND DYNAMICAL 
EVOLUTION 



As stressed by Kandrup and co-workers (Kandrup fc Ma- 



hon 1994, Kandrup 199C), a useful way to think about the 
consequences of chaotic motion for the evolution of stellar 
systems is in terms of mixing. Mixing is the process by which 
an initially non-uniform distribution of points in phase space 
r elaxes to a uniform one, at lea st in a coarse-grained sense 
(Lichtenberg & Lieberman 1992). It therefore provides a link 
between the behavior of individual orbits and the evolution 
of the phase-space density. While there are man y discussions 
of mi xing in the stellar d ynamics litera ture ( Lynden-Bell 
1967| , [Premaine et al. 1994 Ma.thur 198'^ , poker 1996D , most 



of these are vague with regard the origin of the mixing or 
its properties. Kandrup et al. (1994) and Merritt & Valluri 
(1996) noted many of the special properties of chaotic mix- 
ing. It is irreversible, in the sense that an infinitely fine tun- 
ing is required to undo its effects. It has an associated time 
scale, the Liapunov time, in the sense that a compact set of 
points will initially separate at a rate determined by the Lia- 
punov exponent. Following the arguments presented above, 
however, we expect mixing to produce significant changes in 
the macroscopic distribution of phase space points only if 
there are no barriers to the diffusion. 

Figure 4 shows examples of chaotic mixing in a triaxial 



with a shallow cusp and a central point mass. Three en- 
sembles of orbits were started at rest on an equipotential 
surface and integrated in tandem for several crossing times. 
The central point had a mass Mh = 0.03 (in units of the 
galaxy mass M), similar to the ma sses of the largest black 



holes observed in nearby galaxies (Kormendy & Richstonc 



1995). The first ensemble (a) was begun on an equipotential 
surface enclosing a mass ~ 3 times that of the "black hole" ; 
for ensembles (b) and (c) these ratios were ~ 7 and ~ 17 re- 
spectively. The number in each frame is the elapsed time in 
units of the local crossing time. All ensembles were initially 
distributed in a tiny square patch as for T = in ensemble 
(a). 

Mixing occurs very rapidly in these ensembles. At the 
lowest energy, (ensemble a), the linear extent of the points 
in configuration space doubles roughly every crossing time 
until T « 4, when the volume defined by the equipotential 
surface appears to be nearly filled. At the highest energy 
(ensemble c), mixing is slower but substantial changes still 
take place in a few crossing times. The final distribution of 
points at this energy still shows some structure, reminiscent 
of a box orbit. 

Mixing in triaxial potentials with smaller central masses 
can be slower, requiring hundreds or even t housands of cross- 



ing t imes to reach a near-invariant state (Merritt & Valluri 



199C ) - even though the Liapunov exponents in these poten- 



tials are comparable (expressed in units of the local orbital 
frequency) to those of the orbits in Fig. 4. Here again, we 
find that Liapunov exponents are poor predictors of the rate 
of evolution. 

The mixing illustrated in Fig. 4 takes place in a phase 
space that is globally chaotic, with almost no regular orbits. 
In triaxial potentials containing a central black hole, a "zone 
of chaos" extends from the nucleus out to a radi us contain- 
ing ~ 10^ ti mes the mass of the central point (Merritt & 
Valluri 1999). Each of the ensembles in Fig. 4 lies within 



this zone. Near the outer edge of the chaotic zone (ensem- 
ble c), the mixing is beginning to become non-random, and 
the distribution of points at the final time step still shows 
hints of structure. At still larger radii, or in potentials with 
smaller black holes, phase space is a complicated mixture 
of regular and chaotic trajectories (cf. Fig. 3) and mixing is 
correspondingly slower. 

Mixing like that of Fig. 4 should produce rapid changes 
in the shape of a galaxy. Such evolution has in fact been 
observed in A''-body simulations of the response of a triaxial 
gala xy to the growth of a centr al black hole. Merritt & Quin- 
lan (Merritt & Quinlan 1998) found that a triaxial galaxy 



evolves to axisymmetry in little more than the local cross- 
ing time at each radius when the black hole mass exceeds 
^ 2.5% of the total galaxy mass. There are two factors that 
can cause mixing at large radii to occur faster in these A''- 
body simulations than in Fig. 4(c). First, in a self-consistent 
potential all orbits are initially well spread out and not con- 
fined to tiny patches as in the mixing experiments. Second, 
mixing timescales at small radii are almost 100 times faster 
(in absolute time units) than at the half mass radius. Once 
the inner region changes shape the orbits further out are no 
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Figure 4. Mixing in a triaxial potential with a central point containing 3% of the total mass. Time is in units of the local crossing time. 
Ensembles of 10* phase points were distributed initially (T = 0) in patches on an equipotential surface with zero velocity. Ensemble (a) 
was begun on a surface enclosing a mass ~ 3 times that of the "black hole" ; for ensembles (b) and (c) these ratios were ~ 7 and ~ 17 
respectively. Mixing occurs very rapidly for these ensembles, with a characteristic time of order the crossing time. 
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longer in equilibrium in the new potential and can r espond 



mor e ra pidly than they would in a fixed potential (Barne; 
199S|)~|rherefore the orbits at larger radii respond collec- 
tively both to the perturbing central black hole and to the 
potential changes at small radii. Smaller black holes induce 
slower evolution, consistent (at least qualitativel y) with the 
less effi cient mixing expected in such potentials (Merritt & 
Valluri 1996). Such experiments - which should be repeated 



with a wider range of initial conditions, including rotating 
models - confirm that chaos can be an important mechanism 
in determining the global structure of galaxies. 



5 CONCLUSIONS 

The orbits of stars in stellar systems are generically chaotic. 
Exponential sensitivity to perturbations implies that ini- 
tially adjacent trajectories will diverge in a fraction of a 
crossing time. However this divergence does not necessarily 
imply a macroscopic evolution of the phase space distribu- 
tion on a comparable time scale, as has been suggested by 
some authors. The reason is that perturbations may grow 
exponentially only for a limited time, saturating when sepa- 
rations are much smaller than the size of the system. This oc- 
curs, for instance, in a potential composed of point masses 
when A'^ is large. In such potentials, orbits can mimic regu- 
lar orbits for many oscillations even though their Liapunov 
exponents are large. Only when perturbations are able to 
grow until their scale is of order the size of the system does 
the chaos have a significant infiuence on the macroscopic 
dynamics. 

The redistribution of stars in phase space due to chaos 
may be simulated by mixing experiments in fixed poten- 
tials. Mixing leads to a near- invariant distribution of points 
within the accessible phase-space region. We find that mix- 
ing in a globally chaotic region of phase space - a region 
where almost all the orbits are stochastic - is very efficient, 
producing substantial changes in the distribution of points 
in just a crossing time. In phase space regions containing a 
mixture of regular and chaotic trajectories, mixing can be 
slower. The objects that impede mixing appear to be the 
resonant tori and their associated families of regular orbits. 
In triaxial potentials containing a central point mass, mixing 
becomes very efficient when the central mass is big enough 
to destroy all of the resonant tori. The result is large-scale 
evolution of the galaxy toward axisymmetric shapes. 
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